Bayesian statistical modeling by Stan&R introduction (the book by baba)

data was downloaded from this website

library(tidyverse)
library(rstan)
library(bayesplot)
path <- "/Users/tomoya/Documents/dataset/bayesian_intro_baba/2-4-1-beer-sales-1.csv"
df <- read_csv(path)
df
# A tibble: 100 x 1
   sales
   <dbl>
 1  87.5
 2 104. 
 3  83.3
 4 132. 
 5 107. 
 6  83.6
 7 110. 
 8 115. 
 9 112. 
10  93.9
# … with 90 more rows
mcmc_result <- stan(file = "var-mean.stan", data = list(sales = df$sales, N = nrow(df)), 
    seed = 1, chains = 4, iter = 2000, warmup = 1000, thin = 1)

SAMPLING FOR MODEL 'var-mean' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 3e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.053426 seconds (Warm-up)
Chain 1:                0.03644 seconds (Sampling)
Chain 1:                0.089866 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'var-mean' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 9e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.03837 seconds (Warm-up)
Chain 2:                0.035696 seconds (Sampling)
Chain 2:                0.074066 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'var-mean' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 8e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.048032 seconds (Warm-up)
Chain 3:                0.031268 seconds (Sampling)
Chain 3:                0.0793 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'var-mean' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 8e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.048851 seconds (Warm-up)
Chain 4:                0.064571 seconds (Sampling)
Chain 4:                0.113422 seconds (Total)
Chain 4: 
mcmc_sample <- rstan::extract(mcmc_result, permuted = FALSE)
mcmc_sample
, , parameters = mu

          chains
iterations   chain:1   chain:2   chain:3   chain:4
      [1,] 100.64870 104.79868 100.68753 101.61185
      [2,] 102.79747 105.97725 102.40376 101.37497
      [3,] 100.64631 101.67638 101.53096 102.84214
      [4,] 100.56696 101.09195 102.40958 103.42633
      [5,] 102.79871 103.24347 104.02485 103.08640
      [6,] 102.69393 103.24347  99.64743 100.91259
      [7,] 101.57675 102.21287  97.71370  98.74514
      [8,] 102.89931 100.23452 103.89675 104.65942
      [9,]  99.62106  99.40103 101.03884 102.73738
     [10,] 103.55490  99.40626 102.43715 102.94000
     [11,] 102.90797 103.55604 102.74112 102.62739
     [12,] 103.86359  99.09045 104.93339 102.48959
     [13,] 101.24143  99.09045 100.14706 104.38543
     [14,]  99.85707 105.49999  99.87357  99.39213
     [15,] 101.37068 105.18203 100.20999 105.90998
     [16,] 103.57256 104.52472 101.36059 103.34531
     [17,] 101.92233 102.00397 105.07658 100.00039
     [18,] 103.01964 101.83571 101.04252 103.92699

 [ reached getOption("max.print") -- omitted 982 row(s) and 2 matrix slice(s) ]
mcmc_sample[1, "chain:3", "mu"]
[1] 100.6875
mcmc_sample[, "chain:1", "mu"]
 [1] 100.64870 102.79747 100.64631 100.56696 102.79871 102.69393 101.57675
 [8] 102.89931  99.62106 103.55490 102.90797 103.86359 101.24143  99.85707
[15] 101.37068 103.57256 101.92233 103.01964 104.38627 101.61651 103.18445
[22] 100.15041 104.19408 101.60059 100.02435 103.37332 104.80943 101.02734
[29]  99.75564 100.35635 100.97986 101.48021 101.62266 101.94145 101.89881
[36] 102.81174 101.82513 102.36507 101.81274 102.49252 101.28915  99.51839
[43]  99.35342  98.70641 102.85533 101.80770 101.55767 102.99445 101.20891
[50] 102.10344 104.19419 104.39294 101.11485  99.07752 100.25881 100.33359
[57] 100.01780 100.70371  99.51196 100.53404 103.08729 102.13608 103.13025
[64] 106.22945  98.28900 104.21365 103.10044 101.40143 103.48623 103.35901
[71] 104.54869 104.36264 103.02141 102.05826 101.08636
 [ reached getOption("max.print") -- omitted 925 entries ]
mcmc_sample[, , "mu"]
          chains
iterations   chain:1   chain:2   chain:3   chain:4
      [1,] 100.64870 104.79868 100.68753 101.61185
      [2,] 102.79747 105.97725 102.40376 101.37497
      [3,] 100.64631 101.67638 101.53096 102.84214
      [4,] 100.56696 101.09195 102.40958 103.42633
      [5,] 102.79871 103.24347 104.02485 103.08640
      [6,] 102.69393 103.24347  99.64743 100.91259
      [7,] 101.57675 102.21287  97.71370  98.74514
      [8,] 102.89931 100.23452 103.89675 104.65942
      [9,]  99.62106  99.40103 101.03884 102.73738
     [10,] 103.55490  99.40626 102.43715 102.94000
     [11,] 102.90797 103.55604 102.74112 102.62739
     [12,] 103.86359  99.09045 104.93339 102.48959
     [13,] 101.24143  99.09045 100.14706 104.38543
     [14,]  99.85707 105.49999  99.87357  99.39213
     [15,] 101.37068 105.18203 100.20999 105.90998
     [16,] 103.57256 104.52472 101.36059 103.34531
     [17,] 101.92233 102.00397 105.07658 100.00039
     [18,] 103.01964 101.83571 101.04252 103.92699
 [ reached getOption("max.print") -- 982 行を無視しました ] 
mu_mcmc_vec <- as.vector(mcmc_sample[, , "mu"])
mu_mcmc_vec
 [1] 100.64870 102.79747 100.64631 100.56696 102.79871 102.69393 101.57675
 [8] 102.89931  99.62106 103.55490 102.90797 103.86359 101.24143  99.85707
[15] 101.37068 103.57256 101.92233 103.01964 104.38627 101.61651 103.18445
[22] 100.15041 104.19408 101.60059 100.02435 103.37332 104.80943 101.02734
[29]  99.75564 100.35635 100.97986 101.48021 101.62266 101.94145 101.89881
[36] 102.81174 101.82513 102.36507 101.81274 102.49252 101.28915  99.51839
[43]  99.35342  98.70641 102.85533 101.80770 101.55767 102.99445 101.20891
[50] 102.10344 104.19419 104.39294 101.11485  99.07752 100.25881 100.33359
[57] 100.01780 100.70371  99.51196 100.53404 103.08729 102.13608 103.13025
[64] 106.22945  98.28900 104.21365 103.10044 101.40143 103.48623 103.35901
[71] 104.54869 104.36264 103.02141 102.05826 101.08636
 [ reached getOption("max.print") -- omitted 3925 entries ]
mean(mu_mcmc_vec)
[1] 102.166
quantile(mu_mcmc_vec, probs = c(0.025, 0.975))
     2.5%     97.5% 
 98.54126 105.75601 
print(mcmc_result, probs = c(0.025, 0.5, 0.975))
Inference for Stan model: var-mean.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd    2.5%     50%   97.5% n_eff Rhat
mu     102.17    0.03 1.82   98.54  102.19  105.76  2889    1
sigma   18.18    0.02 1.30   15.90   18.09   20.98  3296    1
lp__  -336.44    0.02 0.97 -339.06 -336.15 -335.47  1877    1

Samples were drawn using NUTS(diag_e) at Thu Jan  2 19:00:40 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
library(ggfortify)
autoplot(ts(mcmc_sample[, , "mu"]), facets = F, ylab = "mu", main = "trace plot") + 
    theme_bw(base_size = 14)

ts(mcmc_sample[, , "mu"])
Time Series:
Start = 1 
End = 1000 
Frequency = 1 
       chain:1   chain:2   chain:3   chain:4
   1 100.64870 104.79868 100.68753 101.61185
   2 102.79747 105.97725 102.40376 101.37497
   3 100.64631 101.67638 101.53096 102.84214
   4 100.56696 101.09195 102.40958 103.42633
   5 102.79871 103.24347 104.02485 103.08640
   6 102.69393 103.24347  99.64743 100.91259
   7 101.57675 102.21287  97.71370  98.74514
   8 102.89931 100.23452 103.89675 104.65942
   9  99.62106  99.40103 101.03884 102.73738
  10 103.55490  99.40626 102.43715 102.94000
  11 102.90797 103.55604 102.74112 102.62739
  12 103.86359  99.09045 104.93339 102.48959
  13 101.24143  99.09045 100.14706 104.38543
  14  99.85707 105.49999  99.87357  99.39213
  15 101.37068 105.18203 100.20999 105.90998
  16 103.57256 104.52472 101.36059 103.34531
  17 101.92233 102.00397 105.07658 100.00039
  18 103.01964 101.83571 101.04252 103.92699
 [ reached getOption("max.print") -- 982 行を無視しました ] 
ggplot(data = mcmc_sample[, , "mu"] %>% as.vector %>% data.frame(sample = .)) + 
    geom_density(aes(x = sample), fill = "skyblue3", color = "skyblue4", alpha = 0.6) + 
    theme_bw(base_size = 14)

mcmc_hist(mcmc_sample, pars = c("mu", "sigma"))

mcmc_dens(mcmc_sample, pars = c("mu", "sigma"))

mcmc_trace(mcmc_sample, pars = c("mu", "sigma"))

mcmc_combo(mcmc_sample, pars = c("mu", "sigma"))

mcmc_intervals(mcmc_sample, pars = c("mu", "sigma"), prob = 0.8, prob_outer = 0.95)

mcmc_areas(mcmc_sample, pars = c("mu", "sigma"), prob = 0.8, prob_outer = 0.95)

mcmc_acf_bar(mcmc_sample, pars = c("mu", "sigma"))

animal_num <- read_csv("/Users/tomoya/Documents/dataset/bayesian_intro_baba/2-5-1-animal-num.csv")
animal_num %>% table
.
 0  1  2  3  4 
64 85 39  8  4 
mcmc_normal <- stan(file = "animal_normal_dist.stan", data = list(animal_num = animal_num$animal_num, 
    N = nrow(animal_num)), seed = 1, chains = 4, iter = 2000, warmup = 1000, 
    thin = 1)

SAMPLING FOR MODEL 'animal_normal_dist' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 1.4e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.04218 seconds (Warm-up)
Chain 1:                0.033597 seconds (Sampling)
Chain 1:                0.075777 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'animal_normal_dist' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 1.2e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.05225 seconds (Warm-up)
Chain 2:                0.033195 seconds (Sampling)
Chain 2:                0.085445 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'animal_normal_dist' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 8e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.038179 seconds (Warm-up)
Chain 3:                0.04208 seconds (Sampling)
Chain 3:                0.080259 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'animal_normal_dist' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 8e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.035994 seconds (Warm-up)
Chain 4:                0.035385 seconds (Sampling)
Chain 4:                0.071379 seconds (Total)
Chain 4: 
mcmc_poisson <- stan(file = "animal_poisson_dist.stan", data = list(animal_num = animal_num$animal_num, 
    N = nrow(animal_num)), seed = 1, chains = 4, iter = 2000, warmup = 1000, 
    thin = 1)

SAMPLING FOR MODEL 'animal_poisson_dist' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 1.3e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.035789 seconds (Warm-up)
Chain 1:                0.034036 seconds (Sampling)
Chain 1:                0.069825 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'animal_poisson_dist' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 9e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.033882 seconds (Warm-up)
Chain 2:                0.028672 seconds (Sampling)
Chain 2:                0.062554 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'animal_poisson_dist' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 6e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.031612 seconds (Warm-up)
Chain 3:                0.027575 seconds (Sampling)
Chain 3:                0.059187 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'animal_poisson_dist' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 8e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.028679 seconds (Warm-up)
Chain 4:                0.026872 seconds (Sampling)
Chain 4:                0.055551 seconds (Total)
Chain 4: 
y_rep_normal <- extract(mcmc_normal)$pred
y_rep_poisson <- extract(mcmc_poisson)$pred
dim(y_rep_normal)
[1] 4000  200
y_rep_normal[1, ]
 [1]  1.146980013  1.854437914  1.021145357  0.578413718  0.258020066
 [6] -0.351310680  1.253304848  0.110552641  0.755742723  1.228747160
[11]  1.344334788  0.541145278  0.780954016  1.949987041  0.043828549
[16]  0.721685533  2.040762643  1.010444509  1.047763208  2.621768586
[21]  1.742483477  1.983092709  2.512677485  1.647432538 -0.827800346
[26]  2.267476210  0.883656061  0.629227036  0.947072106  3.769873426
[31] -0.003849530  0.688832818  1.612792468  1.174575527  1.192150488
[36]  2.077392144  1.018682377  0.629081869  2.157806212  1.411999471
[41]  1.154246517  1.852153406  0.779993220  1.153493484 -0.686172286
[46]  1.687172205  1.037327786  0.029628722  1.555700961  1.719985511
[51]  1.215092548 -0.002299914  0.521871637  1.369611267  0.558290493
[56] -0.855950558  1.335117234  2.746824217  0.752569171  1.938582522
[61]  0.588832991  0.801822983  3.518490204 -0.113727357  1.036300359
[66]  0.718570796  0.640913945 -0.140313360  0.621811302  2.642044739
[71]  0.823595780  0.658087706  0.454968067  0.845824208  1.100157136
 [ reached getOption("max.print") -- omitted 125 entries ]
y_rep_poisson[1, ]
 [1] 1 1 2 1 0 0 1 0 1 2 3 0 1 3 1 1 1 0 0 1 1 2 3 6 1 0 1 5 0 0 0 3 0 3 0
[36] 0 4 0 1 1 0 0 2 2 1 1 0 1 1 1 1 0 2 0 0 0 1 0 1 0 1 1 0 0 0 0 3 1 1 0
[71] 2 1 0 1 1
 [ reached getOption("max.print") -- omitted 125 entries ]
ppc_hist(y = animal_num$animal_num, yrep = y_rep_normal[1:5, ])

path <- "/Users/tomoya/Documents/dataset/bayesian_intro_baba/2-6-1-beer-sales-ab.csv"
sales_ab <- read_csv(path)
ggplot(data = sales_ab, aes(x = sales, y = ..density.., color = beer_name, fill = beer_name)) + 
    geom_histogram(alpha = 0.5, position = "identity") + geom_density(alpha = 0.5, 
    size = 0) + theme_bw(base_size = 14)

inspect <- sales_ab %>% mosaic::inspect()
inspect$categorical
# A tibble: 1 x 6
  name     class    levels     n missing distribution                      
  <chr>    <chr>     <int> <int>   <int> <chr>                             
1 beer_na… charact…      2   200       0 "A (50%), B (50%)                …
inspect$quantitative
# A tibble: 1 x 11
  name  class     min    Q1 median    Q3   max  mean    sd     n missing
  <chr> <chr>   <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <int>   <int>
1 sales numeric  55.7  103.   128.  165.  239.  136.  41.1   200       0
data_list_ab <- list(sales_a = sales_ab %>% filter(beer_name == "A") %>% select(sales) %>% 
    1[[]], sales_b = sales_ab %>% filter(beer_name == "B") %>% select(sales) %>% 
    1[[]], N = 100)
mcmc_result_sales_ab <- stan(file = "sales_ab_normal.stan", data = data_list_ab, 
    seed = 1, chains = 4, iter = 2000, warmup = 1000, thin = 1)

SAMPLING FOR MODEL 'sales_ab_normal' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 1e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.050558 seconds (Warm-up)
Chain 1:                0.019665 seconds (Sampling)
Chain 1:                0.070223 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'sales_ab_normal' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 5e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.040216 seconds (Warm-up)
Chain 2:                0.021298 seconds (Sampling)
Chain 2:                0.061514 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'sales_ab_normal' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 5e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.037421 seconds (Warm-up)
Chain 3:                0.04389 seconds (Sampling)
Chain 3:                0.081311 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'sales_ab_normal' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 4e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.047191 seconds (Warm-up)
Chain 4:                0.04551 seconds (Sampling)
Chain 4:                0.092701 seconds (Total)
Chain 4: 
mcmc_sample <- extract(mcmc_result_sales_ab, permuted = F)
mcmc_dens(mcmc_sample, pars = "diff") + theme_bw(base_size = 14)

Tomoya Saito

2020-01-02